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Abstract. We present a phenomenological time-dependent Ginzburg-Landau theory 
of nonlinear plastic deformations in solids. Because the problem is very complex, we 
first give models in one and two dimensions without vacancies and interstitials, where 
large strains produce densely distributed slips but the mass density deviations remain 
small except near the tips of slips. Next we set up a two-dimensional model including a 
vacancy field (or local free- volume fraction) , where relevant is the sensitive dependence 
of the elastic shear modulus on the vacancy density. In our simulation, if strains are 
applied to nearly defectless solids but in the presence of such elastic inhomogeneity, 
the vacancy density and the mass density can become considerably heterogeneous for 
large strains on spatial scales much longer than the atomic size. These strain-induced 
disordered states are metastable or long-lived once they are created. 



1. Introduction 

Most previous papers on glass transitions so far are concerned with near-equilibrium 
properties such as relaxations of the density time correlation functions or dielectric 
response. However, these quantities are too restricted or indirect, and there remains a 
rich group of unexplored problems in far-from-equilibrium states. In particular, shear is 
a relevant perturbation drastically changing the glassy dynamics when the shear rate 7 
exceeds the inverse of the structural relaxation time r a [1, 2]. In supercooled liquids, the 
microscopic rearrangement processes occur on the time scale of r a in quiescent states, 
whereas they are much accelerated even by extremely small 7 (if larger than r^ 1 ). 
Similar jamming rheology has been studied in systems composed of large elements such 
as colloidal suspensions, dense microemulsions, and granular materials [3]. 

Experimentally, Simmons et al. [4] observed shear-thinning behavior, where the 
steady-state viscosity 77(7) was represented by 

»7(T) = ^/7 = '?(0)/(l + 7^) J (1) 

in the range 7x 10 5 < 7/(0) < 6 x 10 13 (Poise) in so da-lime- silica glass. The characteristic 
time r v is expected to be of order r a . Remarkably, the shear stress a xy tends to a limiting 
shear stress, o k™ = r](0)/r v , of order 10 _2 /io, (jlq being the shear modulus for infinitesimal 
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strains. After application of shear, they also observed an overshoot of the shear stress 
before approach to a steady state. 

As a closely related problem, understanding of mechanical properties of amorphous 
metals has been of great technological importance [5, 6]. They are usually ductile in 
spite of their high strength. At low temperatures T <^ 0.6 ~ 0.7T g , localized shear 
bands, where zonal slip occurs, have been observed above a yield stress. At relatively 
high temperatures T > 0.6 ~ 0.7T g , on the other hand, shear deformations are induced 
homogeneously (on macroscopic scales) throughout samples, giving rise to viscous flow 
with strong shear thinning behavior. In a model amorphous metal in 3D, Takeuchi et 
al. [7] numerically followed atomic motions after application of a small shear strain 
to observe heterogeneities among poorly and closely packed regions. Such dynamic 
heterogeneities have been reported in recent simulations in sheared states [1] and also 
in quiescent states [8, 9, 10, 11]. 

In the study of elasto-plastic dynamics of solids, microscopic simulations are 
informative [12, 13, 14], while mean-field theories are instructive [15, 16]. In the latter 
theories, the problem was reduced to to that of one element obeying a stochastic process 
under the influence of the average stress. The aim of this paper is then to present a space- 
time dependent, elasto-plastic dynamical model on the basis of well-defined nonlinear 
elastic theory. Although our theory is still preliminary, we shall see that introduction 
of an order parameter, representing the vacancy or the local free volume, can give rise 
to dramatic effects similar to those reported in Ref.[12]. 

2. Plastic flow in one dimension 

To introduce the fundamental concepts in our problem in the simplest manner, we first 
present a one-dimensional model which mimics a solid with shear deformations varying 
only in one direction (ID slip model). We here write down the dynamic equations in 
the continuum representation (for simplicity). The velocity v(x,t) = du(x,t)/dt of the 
shear displacement u(x,t) (along the y axis) is governed by 

d 

P dt V = ^ x ^° 7 ° sin( ^ 7o ~ 1 ^ xU ^ + Vo^lv + V x a R , (2) 

where V x = d/dx, p is the mass density, fio is the shear modulus, 70 is a constant 
determined by the underlying crystal structure, and r/o is a vicosity. The period of the 
strain is given by 7 P = 27T7 . To ensure equilibrium in the absence of applied force, 
we introduce the random stress <tr(x, t) related to i]q via the fluctuation-disipation 
relation (or(x, t)a R (x', t')) = 2/c B T?7o<5(a; — x')5(t — t), where T is the temperature. 
Hereafter space and time will be measured in units of i = ^/(p/Uo) 1 / 2 and Uq 1 = i] //j,q, 
respectively. The dimensionless noise strength is given by e = k^I '/ '/io7o^- 



3 



We numerically solve (2) on a ID chain (j = 0, 1, • • • , N = 600) with e = 0.25, 
setting u(x,t) — > uj(t) and du/dx — > jj = Uj +1 — Uj. We apply a constant shear rate 
7 at t = with uq = and ttjv/7o = Njt. In Fig.l the average dimensionless stress 
(<j)(t) = iV _1 J2j sm (7j(^)/7o) + 7 is plotted as a function of the scaled strain A ft/ r y p . 
The system undergoes a pure elastic deformation 7j(£)/7o — 7^ in the initial stage, 
but slips (jumps of 7j by multiples of 7 P ) appear randomly throughout the system with 
increasing the average strain jt. The inset shows slips and continuously strained regions. 
The average slope of the latter regions gives the average stress. For the range of the 
shear rates in Fig.l, we can see shear-thinning a xy oc 7°' 7 in steady states (if averaged 
over time). 

It is also straighforward to integrate (2) under a fixed shear stress <7 ex t applied at 
one end with the other end being pinned. For small <7 cx t below a yield stress a y , we 
observe no slip formation, while for <7 ext > a y plastic flow is produced. For the noise 
strength of e = 0.25, we find a y = 0.2. 



3. Plastic flow without vacancies 



As a direct generalization of (2), we set up a plastic flow model in two dimensions (2D 
slip odel). In the continuum limit the strain components are defined by 

e 1 = V x u x + VyUy, e 2 = V x u x - V y u y , e 3 = V ' x u y + V y u x , (3) 

where V z = d/dx and V y = d/dy. If we suppose a triangular lattice, the elasticity is 
isotropic in the harmonic approximation [2], being characterized by the bulk and shear 
moduli, Kq and fiQ, but it depends on the orientational angle 6 of one of the crystal axes 
for large shear deformations. Note that, under rotation of the reference frame by 9, the 
shear strains e 2 and e 3 are changed to e' 2 and e 3 , where [2] 

e' 2 = e 2 cos 29 + e 3 sin 29, e' 3 = — e 2 sin 29 + e 3 cos 29. (4) 

We write the elastic energy density in the form, f e \ = K el/2+fi J : '(e' 3 , e' 2 ). The simplest 
form of the scaling function JF is of the form 



jF(e 3 ,e 2 ) = -^-^ 3 — cos7r(v / 3e 3 — e' 2 ) — cos7r(v / 3e 3 + e' 2 ) — cos(27re 2 ) 



(5) 



and is shown in Fig. 2, The T is invariant with respect to the rotation 9 — > 9 + 7r/3, is 
a periodic function of e' 3 with period 7 P = 2/a/3 for e' 2 = (simple shear deformation), 
and becomes (e 2 + e 3 )/2 for small strains. If we assume 89/ dt = (V x v y — V y v x )/2 for 
the angle rotation rate, its time integration yields 

= \iyxU y -V y u x ) + 9 , (6) 
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where 9 is independent of t but may depend on r = (x, y). For 9 = one of the crystal 
axes is along the x axis. Then the elastic energy F el = / drf c i is a functional of u. We 
assume that the lattice velocity v = du/dt obeys 

d & 

p-^v = -—F el + r) V 2 v + V-a R , (7) 

where the first term on the right hand side is also written as V- o~ in term of the elastic 
stress tensor a. The symmetric random stress tensor <j r = {cr^p} = {°~p a } satisfies 
af x + a^y = 0, because the bulk viscosity is neglected in (7), and [2] 

(a^(r,t)a^(r',t')) = 2k B Tr) 6(r - r')5(t - t'). (8) 

We measure space and time in units of ^/(p/A)) 1 ^ 2 and i] / p an d the strains in 
units of 7 P = 2/y/3. If these scaling units are used, the noise strength k B Tr] in (8) is 
replaced by 

e = k B Tphlr,l. (9) 

We integrate (7) on a 128 x 128 square lattice by applying a constant shear rate 7 at 
t — with e = 0.1. The periodic boundary condition is imposed in the x direction, 
while u = at the bottom y = and u x = 'yt and u y = at the top y = 128. At t — 0, 
the values of v at the lattice sites are Gaussian random numbers with variance e 1//2 but 
u = and 6 = 0.1. In Fig. 3 we show the average scaled stress (cr xy ) / fi as a function 
of the scaled strain 7t/7 p for 7 = 10~ 3 and 10~ 4 . In the inset we display snapshots of 
e3 and e<i- Because we start with a perfect crystal with the initial fluctuations only in 
the velocity, the stress curve drops sharply after the peak with catastrophic formation 
of slips. Then structurally disordered states are produced where defects are proliferated 
(strain-induced disordering). Note that slips are accumulations of dislocations [17]. In 
the initial stage of plastic flow we have simple slips consisting of two dislocations with 
opposite Burgers vectors ±a with a being the lattice constant. The elastic energy to 
create such a slip is minimum in the x and y directions under shear deformation. This is 
the reason why the slips in Fig. 3 are parallel to the x or y direction. In Fig.4 the shear 
(= 10~ 3 ) is swiched off at (a) t = 190 before the peak time of the stress, (b) t = 220 just 
after the peak time, and (c) t = 440. The top and bottom boundaries are kept at rest 
afterwards. Pure elastic deformation is maintained in (a), while no appreciable time 
evolution is detected after transients in (b) and (c). This means that the structurely 
disordered states are metastable. 

4. Plastic flow with vacancies 

Usually in the literature, the density deviation 5p is equated with — pV ■ u = —pe\ 
in solids [17], where p(= p) is the average density and \5p\ <C p is assumeed. In the 
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presence of vacancies or interstitials, however, there can be a small difference between 
these two quantities. Cohen, Flemming, and Gibbs [18] constructed a hydrodynamic 
description of solids including the new variable m = Sp/p + e±. According to their 
theory, the vacancy concentration c may be defined by 

c = c - m = c - (Sp/p + ei) (10) 

where Co is the average (c) dependent on T and p. On the other hand, Granato claimed 
relevance of interstitials in amorphous solids because they can greatly decrease the shear 
modulus [20]. However, because vacancies and interstitials are point defects, it may be 
more appropriate to treat c as the local free-volume fraction [18, 19] which can take 
continuous values at each lattice site. Then Co is the average free-volume fraction and 
we may assume < c C 1. The role of m or c is expected to be crucial in amorphous 
solids. We try to include the vacancy variable in our nonlinear slip model (2D vacancy 
model), 

The dynamic equations obeyed by the mass density and the momentum density are 
of the usual forms, 

sLs P = -pv-v, (ii) 

d 

P^2 V = V- ? + Vo V 2 v + W-a R . (12) 
ot 

In (12) the force density, the first term on the right hand side, is written in terms of the 
free energy functional F = F{p, u} as 

V- a= -pV(5F/5p) u - (SF/5u) p = -{5F/5u) m , (13) 

where use has been made of the identities (5/8u) p = (5/5u) m — V(5/5m) u and 
(5/8m)u = p{S/Sp)u- The random stress tensor <x R satisfies (8). It is important 
that, in the presence of the vacancy field, the lattice velocity du/dt is different from the 
mass velocity v and is assumed to be of the form, 

| u = t ,_ Ao (_L F ) /<R , (14) 



where Ao is the kinetic coefficient and the components (^ °f the random force vector £ R 
are characterized by 

(O, t)(f(r', 0) = 2k B T\ 5 a p5(r - r')5(t - t'). (15) 

From (12) and (14) the equation for m is expressed as 

+ V-Cr- (16) 



9 ^ \ 
—m = V • A 
ot 



v(—F) - Wt 

\om 



u 
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It is worth noting that (16) is similar to the dynamic equation for the concentration in 
visccoelastic fluid mixtures [2, 21]. In passing, owing to the special form of (13), the time 
derivative of the total free energy F{p, u} + / drpv 2 /2 becomes nonnegative-definite in 
the absence of applied stress if the random noises are neglected. This is a self-consistent 
condition of Langevin equations ensuring attainment of equilibrium [2]. 

The free energy F — J drf is a functional of u and Sp (or m). We assume the free 
energy density / in the form, 

/ = ^rn 2 + |m 4 + ame, + +j\V5p\ 2 + ^-e\ + p(m)F(e' 3 , e' 2 ), (17) 

where A, B, and C are positive constants and a is a coupling constant. The gradient 
term (oc C) is introduced to suppress the density fluctuations with short length scales. 
(We are considering the density fluctuations longer than the peak distance of the pair 
correlation function.) For simplicity, we neglect the other gradient terms involving the 
gradients of m and the strains. In our theory it is most important that the shear modulus 
p(m) sensitively depends on m. Near glass transitions, with decreasing m or increasing 
c, pirn) is expected to decrease abruptly from a finite value po to zero (= fluid value) 
around a threshold value of m or c. 

Analytic calculations of the above model are difficult except for idealized situations. 
As a simple exercise, let us consider small deviations with wave vector k around a 
homogeneously strained state with e 3 = 7, where the sytem is at rest and the harmonic 
approximation can be made for the elastic energy. Then the mechanical equilibrium 
condition allows elimination of u from (16) [2]. We consider the steady-state variance of 
m in the long wavelength limit, Xv = hm*;^o(| m fo| 2 ), m ^D. In the linear approximation 
it becomes dependent on the angle of k: 

— = A - ^ - % l7 sin2y? - ^(pa) 2 [l + (cos2v?) 2 7V/io] ; ( 18 ) 

Xv Lq Lq 

where sin2(/? = 2k x k y /k 2 , p(m) = p + p\m + • • • and Lq = K Q + p . The vacancy 
diffusion constant is written as D v = \Qk\ rmB T/xv The steady-state density variance 
X P = linifc^o(|Pfc| 2 )/p 2 m 2D is written as 



k B T , 
X P = -j— + 
Lo 



1 + -!- (a + A*i7 sin 2tp) 



1 2 



Xv (19) 



The first term is the usual term, while the second term arises from the vacancy 
fluctuations. The dependence of these quantities on 7 and the angle tp can become 
appreciable with increasing p\ and/or decreasing A. For a — 0, for example, the 
homogeneous state is unstable for 7 > 7c = (pqA) 1 ^ 2 / p 1 against the fluctuations with 
k along the x or y axis. Note that ) c < 1 can hold if pi/po ^> 1. It is easy to 
confirm that essentially the same result follows in 3D in both shear and elongational 
deformations. Thus, in elastically deformed amorphous solids (before onset of plastic 
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flow), we predict enhancement of the density fluctuations in particular directions of the 
wave vector. 

We then show our numerical results of our model in (11)-(17) obtained on a 128 x 128 
square lattice for the case of applying a constant 7. Again the units of space, time, and 
strains are ?7o/(p/io) 1//2 , Vo/f^o, and 7 P = 2/^/3, respectively. The to is also divided by 
7 P . The dimensionless mass density is given by 5p/p7 P . Though it is not clear what 
functional form of //(to) is appropriate, we tentatively use 

H(m) = fi exp[-A /(m + m )] (20) 

for m + m > and //(to) = for m + m < 0. Here A is a constant and is set equal 
to 0.1 in our simulation, while to is a control parameter representing the closeness to 
the glass transition. In fact, if //(to) goes to zero at c = c c , we have m = c c — c from 
(10). However, essentially the same numerical results follow as long as //(to) strongly 
depends on to (for example, from the simpler form G{m) oc (to + to ) 2 instead of (20)). 
The other parameter values used are A — C — 1, B — 10, a — 0, Ao = 0.1 and 
e = 0.1. While the initial conditions for v, u and 9 are the same as in the previous 
case of Figs. 3 and 4, we assign random Gaussian randum numbers with variance 0.1 
to to. In Fig. 5 we show time evolution of the scaled shear stress as a function of the 
scaled average strain 7t/7 P at 7 = 10~ 3 for m = 1,0.6, and 0.2. Fig. 6a and Fig.6b 
display snapshots of e^, e2, c — c = —to, and 5p at t = 250 and 420 on the curve 
of m = 0.6 in Fig.5. Fig. 7 demonstrates that these strain-induced disordered states 
are metastable as in Fig. 4. Salient features are as follows: (i) Plastic flow is induced 
earlier than in the absence of vacancies, because slip formation is easier in regions with 
smaller m. Onset of plastic flow is sensitive to the initial randomness of m. (ii) The 
effective viscosity ?7efr = (cfxy)/i becomes smaller with decreasing the parameter mo in 
(17) or increasing the elastic inhomgeneity. (iii) The 77^ also decreases as a function of 
time on long time scales because of slow accumulation of vacancies around slips. In our 
model, though this tendency is considerably suppressed by the quartic term (oc B) in 
(16), vacancy accumulation can lead to fluidization at long times after percolation of the 
slips, (iii) The density fluctuations are induced and are frozen after cessation of shear. 
As demonstared in Fig. 7, the fluctuation variances of the vacancy and the density are 
much enhanced than in the initial state. Owing to the gradient term in (16), the density 
smoothly varies in space (compared with to). 

5. Concluding remarks 

We have presented a first time-dependent Ginzburg-Landau theory accounting for 
inhomogeneous nonlinear elastic deformations. It is a coarse-grained theory and the 
sharp peak of the structure factor, which is essential in the mode coupling theory [22], 
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does not come into play. Instead, periodicity of the elastic energy with respect to 
the strains e' 3 and e' 2 and vacancy-dependence of the shear modulus are two major 
ingredients of our theory, which bring about metastable structural disordered states 
with heterogeneities in the vacancy and mass densities. As a closely related effect, 
Fischer [23] observed excess scattering from nearly static density fluctuations with sizes 
in a range of 20-200 nm in glass-forming fluids. We believe that such large-scale frozen 
fluctuations can arise only as a result of nonlinear elastic deformations induced by 
structural disorder. We also mention previous simulations. Argon et al. [12] applied 
a tensile (elongational) strain to a 2D model to produce slips and shear bands making 
angles of ±45° with respect to the stretched direction in agreement with experiments. 
The same results also follow from our model and will be reported shortly. We will show 
that the elastic energy to create a slip under tension is minimum for these directions. 
Ikeda et al. [14] applied a tensile strain to a 3D model to induce a change from a perfect 
crystal to an amorphous solid. 

However, we admit that our results are still preliminary and it remains unclear 
how our theory corresponds to real physical effects. In particular, the initial conditions 
(it = and 9 =const.) should be inappropriate for amorphous solids (where structural 
disorder should preexist even before applying strain). In future work, we should start 
with initial states with various amounts of disorder and should also examine what 
parameter values are appropriate for amorphous solids. 

We mention analogous and instructive examples. In solids undergoing phase 
separation or structural transitions, the dependence of the shear modulus on the order 
parameter, written as m, in the form /j, — hq + pu\m is of crucial importance in phase 
ordering [2, 24]. In all these cases, local minimization of the nonlinear elastic energy 
{{5F/8u) m = 0) gives rise to heterogeneous metastable states where time-evolution is 
pinned. In viscoelastic polymer systems, the composition-dependence of the viscoelastic 
properties can give rise to shear-induced phase separation [2, 21], where the scattered 
light intensity takes a form similar to Xv in (18). 
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Figure Captions 

Fig.l, Stress vs strain in units of /i 7o and 7 P obtained from the discretized version of 
the ID model (1) after application of constant shear rate. Here 7/cJo = 10~ 3 , 10™ 4 , 10 -5 , 
and 1CT 6 from above. The inset displays a snapshot of the deviations Wj/70 — jtj 
(0 < j < 600) at t = 720 for 7 = 10~ 3 . 

Fig. 2, Function T[e-i^e2) defined by (5), which represents the 2D shear deformation 
energy density divided by /io for 9 = 0. 

Fig. 3, Stress vs strain for 7 = 10~ 3 and 10~ 4 obtained from the 2D slip model (7). In 
the inset e 3 and e 2 at t — 400 for 7 = 10~ 3 are shown. 

Fig. 4, Scaled stress and elastic energy vs time, which demonstates freezing of strain- 
induced diordered states of (7) at long times. The shear is switched off at (a) t = 190, 
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(b) t = 220, and (c) t = 440, as indicated by arrows, using the same parameters as for 
the curve of 7 = 10~ 3 in Fig. 3. 

Fig. 5, Stress vs strain for 7 = 10~ 3 obtained from integration of (11), (12), and (14) 
using (16) and (17) at 7 = 10~ 3 . Here m = 1,0.6, and 0.2 from above. 

Fig. 6, Snapshots of e 3 , e 2 , c — c = —m, and 5p at (a) t = 250 and (b) 420 for 7 = 10~ 3 
on the curve of m = 0.6 in Fig. 5. 

Fig. 7, Freezing of strain-induced diordered states of the vacancy model. Shown are the 
density variance [((^p) 2 )] 1 ^ 2 /p7 P , the vacancy variance [{(5m) 2 )] 1 / 2 /^ p , and the average 
elastic energy (/ e i). As in Fig. 4, the shear rate is switched off at t = 220. The frozen 
density profile at t — 2 x 10 4 is shown in the inset. 
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